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Abstract 

We study the relationship between local and global error in Runge- 
Kutta methods for initial-value problems in ordinary differential equa- 
tions. We show that local error control by means of local extrapolation 
does not equate to global error control. Our analysis shows that the global 
error of the higher-order solution is propagated under iteration, and this 
can cause an uncontrolled increase in the global error of the lower-order 
solution. We find conditions under which global error control occurs dur- 
ing the initial stages of the RK integration, but even in such a case the 
global error is likely to eventually exceed the user-defined tolerance. 
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1 Introduction 

Runge-Kutta (RK) methods are the most popular choice of one-step methods 
for solving problems of the form 

y' = f(x,v) (i) 

y(xo) = 2/o 

numerically. In the implementation of these methods, local error control via 
local extrapolation is the preferred choice of error control. It is known, however, 
that this form of error control does not amount to control of the global error. In 
this paper, we seek to give our own interpretation and discussion of this problem, 
including an analysis showing that global error control, if it does occur, is limited 
under RK iteration. 
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2 Relevant concepts 



We now define local and global errors in RK methods formally, and study the 
propagation of local error in the implementation of an explicit RK method. We 
will also find the relationship between local and global error. 

2.1 Runge-Kutta methods 

The most general definition of a Runge-Kutta method pQ is 

k p = f I x t + c p h, Wi + hJ2 a pq k q ) p = 1,2, ...,m 

V « =1 / (2) 

w i+ i = Wi + hJ2 bpk p 
p=l 

Such a method is said to have m stages (the k q ). We note that if a pq = 
for all p ^ q, then the method is said to be explicit; otherwise, it is known as 
an implicit RK method. The number of stages is related to the order of the 
method. The symbol w is used here and throughout to indicate the approximate 
numerical solution, whereas the symbol y will denote the true solution. 

We denote an explicit Runge-Kutta method of order z (RKz) for solving (JTJ) 

by 

<+i =w*+hF z (xi,wf) 

where wf denotes the numerical approximation to y (xi) and F z (x, y) is a func- 
tion associated with the particular RKz method. Indeed, F z is simply the linear 
combination of the stages for that particualr method, as in 

m 

F z =^b p k p . 
P =i 

2.2 Local and global errors 

We define the global error in a numerical solution at x i+1 by 



Af +1 = < +1 -y i+ i, 
and the local error at x^+i by 

e z i+l = [yi + hF z (xi, yi )}-y i+1 . (3) 

In the above, yt denotes the true solution y{xi) , and similarly for yi+\. Note 
the use of the exact value yi in the bracketed term in ©. 
Previously, we have shown [2] that 

Af +1 = e! +1 +«fA? (4) 
of = l + hF'ix^i-h), 
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where € {ViiVi + Af). Equation (U]) provides the relationship between local 
and global errors in RKz. We will assume that Ao = (i.e. the initial value is 
known exactly). We see that the global error at any node Xi + \ is the sum of a 
local error term and a term incorporating the global error at the previous node. 
For convenience, we will drop the argument from F z (xi,^; h) in the remainder 
of the paper; its presence is implied. 
For RKz, it is well-known that 

4+i oc h z+1 
A z +1 ex h*. 



2.3 Local error estimation 

Consider two RK methods of order z and z+l. Let w z +1 denote the approximate 
solution at Xi+i obtained with the order z method, and similarly for w *i . Let 
the local error at Xi + \ in the order z method be denoted by e z +1 — (3 z +1 h z+1 , 
and similarly for e z ^_l = f3 z ^h z+2 . Hence, with w\ ,w z+1 = yi, we have 

if h is sufficiently small. This gives 

Pi+i J^+l • ^' 

2.4 Local error control 

Once we have estimated the local error, we can perform error control. Assume 
that we require that the local error at each step must be less than a user-defined 
tolerance 5. Moreover, assume that, using stepsize h, we find 



-i+i | 



\P! + ih z+1 \ >S 



In other words, the magnitude of the local error e z +1 exceeds the desired toler- 
ance. We remedy the situation by determining a new stepsize h* from 



ft +1 [h*]' +1 = S ^ h * = {]j^\) (6) 

and we repeat the RK computation with this new stepsize. This, of course, 
gives 

Xi+i — Xi + h* . 

This procedure is then carried out on the next step, and so on. Such form of 
error control is known as absolute error control. If the estimated error does not 
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exceed the tolerance, then no stepsize adjustment is necessary, and we proceed 
to the next step. 

Often, we introduce a so-called 'safety factor' cr, as in 




where a < 1, so that the new stepsize is slightly smaller than that given by 
^j. This is an attempt to cater for the possibility that f3 z i+ i may have been 
underestimated, due to the assumptions made in deriving ([5]) . The choice of the 
value of cr is subjective, although a representative value is 0.8. 

Note that because this error control algorithm is applied on each step, we 
could find that over the interval of integration we have stepsizes of varying 
lengths. For this reason, it is appropriate to make the replacement 

h^hi = Xi+i - Xi 

in ©. 



2.5 Propagation of the higher-order solution 

There is a very important point that must be discussed. The method for deter- 
mining (3 Z +1 hinged on the requirement wf,w% + — yi. However, we only have 
the exact solution at the initial point xo; at all subsequent nodes, the solution 
is approximate. How do we meet the requirement wf,w z+1 = yil 

In the case of local extrapolation, the answer is simple: simply use the higher- 
order solution w z+1 as input to generate both wf +1 (with the order z method), 
and wl^l (with the order 2 + 1 method). In other words, we are assuming 
that w z+1 is accurate enough, relative to wf , to be regarded as the exact value, 
an assumption entirely consistent with the assumption made in deriving ([5]). 
This means that we determine the higher-order solution at each node, and this 
solution is used as input for both methods in computing solutions at the next 
node. The question of whether or not the global error that accumulates in the 
higher-order solution affects the calculation of j3 z +1 in (JSJ is addressed in the 
next section. 



3 Analysis 
3.1 The problem 

Now, as per the last paragraph of the preceding section, assume that w z+1 is 
used to generate wf +1 and . Such value of wf +1 (and associated quantities) 
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z,z+l 

y l+ l • 



Hence, we have 



2,2+1 _ nZ UZ+1 I 2,2 + 1 A Z+l 



A2.2-I-J- _ QZ , 

^i+l — Pi+l n 

A 2+1 _ OZ+1.2 + 2 , 2+1 A 2 + 1 

— Pt+1 ^ + a t A 



Thus, 

<fi +1 -<i = /9f + i^ +1 + A i * +1 + fc i ^* +1 Af+ 1 

- (#+i >* +2 + A^ +1 + hF* +1 A* + ') 
= PUihi +1 -PUlK +2 + (F:' z+1 -F* +1 )h i A* +1 (7) 

for small ^, because /i*A? +1 = (/i^+ 2 ) . We see that the presence of global 
error in the higher-order solution does not affect the expression for I3 z i+1 obtained 
under the assumption wf,u>l +l = yi, particularly if hi is small. 

However, the expression for A*^ 1 informs of a potential problem: we have 

A^+ 1 =/3f +1 / l r 1 + «r + W +1 , (8) 

where A^ +1 is the global error in w* +1 . In ([7]), we see that a subtractive cancel- 
lation ensures that the A^ +1 term does not enter directly into the estimate for 

/3f + i- Nevertheless, even if \/3* +1 h* +1 1 ^ S, we could still have A^f* 1 > 5, per- 
haps substantially so, if |A^ +1 | is large. Moreover, we should certainly expect 
that |A^ +1 | could become large under iteration (i.e. as i increases), since global 
error is essentially an accumulation of local errors. The point here is that, even 
if local error control is effective, the global error A*+^ could become large, and 
could grow in an uncontrolled fashion. 

3.2 Bounded global error via local error control? 

Let us investigate the effect on the global error if local error control, via local 
extrapolation, is implemented. Consider the expression obtained previously for 
the global error at Xi+i 

Af +1=£ f +1 + afAf. (9) 
We assume Ao = 0. If we have the exact value yi at each node, then we have 

— £ i+l 

at each node, so that the global error is equal to the local error. If the local 
error has been controlled (subject to tolerance §), we have 

Af +1 *S S 
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which means that the global error satisfies the tolerance S. 

However, as discussed previously, we do not have yi at each node. Rather, in 
the case of local extrapolation, where we have a higher-order solution available, 
and we propagate this higher-order solution, we have, from ([5]). 



A 2,2+1 
^l + l 



z+l 



-z+l 



; 1 + hR 



z,z+l\ 



.z+l 



z, z+l z+l 



£ m+E £ f 1 +°(^ +3 ) 



2,2 + 1 2,2 + 1 



2,2+1 Z + l 



< 



H+l 



E^ +2 

i=i 



- 1 hZ+2 



L /3 z +1 h z+1 + 



z+2 



if 



x ) h 



z+l 



where h is assumed to be the stepsize determined from ^ and we have included 
the safety factor a explicitly. For ease of analysis, we assume here a uniform 



stepsize h. In the second last line, we assume that, since e z +1 ^ S and e z+1 <C e 



(the fundamental assumption in local extrapolation), we must have e 



z+l 



— z + l 

We denote the average value of (3 Z +1 on [xq, Xi) by (3 , and we use ih = 



^ 5. 



-x - 



Now, assuming 



f3 z+1 {xi 



A 2,2 + 1 
^l + l 



xq) h z+1 < S (see the Appendix), we have 



w a z+1 S- 


-a z+2 


^ a z+1 5- 


-a z+2 






= (a z+ U 


_ a z+2 



(3 Z+1 (Xi - Xq) 



It is easily confirmed that for a 
This means that 



= 0.8, we have (a z + 1 + a z+2 ) 
A^ 1 - <>■ 



(10) 

< 1 for z > 2. 



for these values of a and z, which suggests that the global error, like the local 
error, satisfies the user-defined tolerance. In other words, propagation of the 
higher-order solution in local error control via local extrapolation, has resulted 
in control of the global error, although the significance of the safety factor in 
deriving this result should be clear. Most importantly, the above result holds 



2,2+1 

H+l 



is 



only if the assumptions made here are true; if they are not, then 
probably greater than S. For this reason, we say that the global error is possibly 
bounded, but this is not guaranteed. We should appreciate that such a bounding 
of the global error is a beneficial by-product of local error control, and is not the 
designated objective. 

The most important assumption made above is 



52+1 



(Xi - Xq) h z+i < 5 
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In the Appendix, we show that this is, in fact, a consequence of the condition 



i+l 



h 



z+l 



> i 



(ii) 



If this condition is violated, then the assumption does not hold. It is worth ex- 
amining this condition in closer detail: the factor i represents iteration number. 
It is quite reasonable to assume that there exists a value of i such that, for the 

- z + l 

given values of p and h, we will have 



< i 



z+l 



,z+2 



In other words, eventually the RK iterative process will cause (ITT1) to be vio- 
lated. This means that sooner or later the global error will exceed the imposed 
tolerance, despite local error control, and that the bounding of the global error 
will be 'short-lived', so to speak. 

We also note that our assumption of a uniform stepsize is reasonable if the 
stepsize does not vary considerably; nevertheless, in reality it may do so, which 
would also compromise the validity of ([TO . 



4 Comments 

Some comments are appropriate: 

1. We are not restricted to using a method of order z + 1 as the higher-order 
method in local extrapolation. Any method of order z + r.r 1 can be 
used. Our analysis and results are essentially unchanged, save for dTUJ) , 
which becomes 

|A*f+ r | < (a z+1 +a z+r+1 )8. 

The value of r here will influence the value of the safety factor a for which 
the coefficient is less than unity. 

2. Our analysis clearly shows that global error control cannot be achieved 
through local extrapolation. Global error control is usually achieved through 
reintegration - estimating the global error after local extrapolation, and 
then redoing the computation on the entire interval of integration with 

a smaller stepsize. This approach, while effective, can be inefficient and 
probably cannot be implemented for real-time problems, where a globally 
accurate result is needed before the next iteration. It is not our inten- 
tion to report on methods which address this issue, but we would like to 
take the opportunity to refer to our own recent work in this regard, where 
we have developed an algorithm based on high-order quenching to enable 
stepwise global error control 3 [4 . 
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5 Numerical example 



An instructive example is the initial- value problem 

, /In 1000 \ 

y(o) = i 

on [0, 100] . The coefficient in the differential equation has been chosen so that 
y does not exceed 1000 on [0,100], i.e. y does not vary substantially, so that 
absolute error control is suitable. The exact solution is y (x) = e '"ioo 00 x . 

We use RK3 [5 and RK4 [5] to implement local error control, with 8 = 10~ 8 . 
The quantities |e 3 | and |a 3 ' 4 A 4 | are shown as functions of x in Figure 1 (figure 
follows appendix). We see that the local error \s 3 \ is bounded by 8, as expected. 
However, |a 3 ' 4 A 4 | increases monotonically. The global error |A 3,4 | is given by 
\e 3 + a 3 ' 4 A 4 1 and, although not shown, it is easy to see that | A 3 ' 4 | must increase 
beyond 8, and is almost 100 times greater than 8 at x = 100. In fact, |A 3 ' 4 | 
becomes larger than 8 at x = 23.6, after 118 iterations. This is the value of 
i for which (fTTj) is violated for this example. For x < 23.6, the condition is 
not violated and the global error is bounded by 8. Clearly, though, this state 
of affairs does not last, and the propagation of A eventually compromises the 
global accuracy of the solution. 

6 Conclusion 

We have investigated the relationship between local error and global error in RK 
methods, under the implementation of local error control via local extrapolation. 
We find that, even though local error is successfully controlled, we cannot expect 
the same for global error. Our analysis shows that there is a possibility that the 
global error will be bounded by a user-imposed tolerance during the initial stage 
of the integration, but this will not last. The propagation of the global error in 
the higher-order method will eventually cause the global error in the lower-order 
method to exceed the tolerance. A numerical example clearly illustrates these 
points. 
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7 Appendix 

7.1 The assumption 



[3 Z+1 (xi - x ) 



h z+l < 5 



Assuming a constant stepsize h, the underlying premise of local extrapolation 
is the assumption that 



# +1 h* +L » Pt+i \h 



2 + 11 h z + 2 



which implies 



f3 z i+1 \h z+1 =M\(3^\h 



2 + 1 I .2 + 2 



where M > 1 is a large number. Assuming that f3*+l is a slowly varying 



function of x, we can replace ftl+l with its average value /3 Z+1 , and so we may 
write 



h z+2 = 



\PUi\h z+1 >i 
subject to the condition that i < M . Hence, 

\(3 z +l \h z+1 < 8 ^ p z+1 { Xi -x ) 



P Z+1 (xi - x ) 



,2 + 1 



< s. 
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Figure 1: Error components \e 3 | and loc 3 ' 4 A 4 l for the example. 

3 4 

Global error IA ' I (not shown) is the sum of these 
two components. 



